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ABSTRACT 


We  present  a  preliminary  model  of  the  three-dimensional  seismic  structure  of  the  Iran  region  obtained  via 
simultaneous,  joint  inversion  of  body  wave  travel  time  and  gravity  observations.  The  body  wave  data  set  is  derived 
from  location  calibration  efforts  and  includes  a  large  (>1000  events)  subset  of  GT5-level  events.  The  arrival  time 
data  sets  for  these  events  include  many  readings  of  direct  crustal  P  and  S  phases,  as  well  as  regional  (Pn  and  Sn)  and 
teleseismic  phases.  These  data  have  been  groomed  to  identify  and  remove  outlier  readings;  empirical  reading  errors 
are  estimated  for  most  arrivals  from  multiple  event  relocation  analysis. 

We  use  both  free-air  and  Bouguer  gravity  anomalies  derived  from  the  global  gravity  model  of  the  GRACE  satellite 
mission.  The  gravity  data  provide  information  on  shallow  density  variations  with  short  spatial  wavelengths,  and 
deeper  density  structures  with  longer  spatial  wavelengths.  To  increase  the  usefulness  of  the  gravity  data,  we  apply 
high-pass  filtering,  yielding  gravity  anomalies  that  possess  higher  resolving  power  for  crustal  and  lithospheric 
structures.  To  integrate  both  data  sets  into  a  single  inversion  a  functional  relationship  between  seismic  velocities  and 
density  is  required.  We  first  use  a  combination  of  two  empirical  relations:  one  most  appropriate  for  sedimentary 
rocks,  and  a  linear  Birch’s  Law  that  is  more  applicable  to  basement  rocks.  These  relationships,  however,  assume  a 
constant  Poisson  ratio  throughout  the  inversion.  We  therefore  also  apply  a  method  that  allows  the  Poisson  ratio  to 
vary  by  mapping  densities  to  VP  and  Vs  independently  for  different  earth  materials. 

Final  results  of  the  simultaneous  inversion  will  help  us  to  better  understand  the  crustal  and  lithospheric  structures 
associated  with  early-stage  continental  collision.  Such  models  also  provide  an  important  starting  model  for 
computationally  more  expensive  and  time-consuming  fully  3D  waveform  inversions. 
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OBJECTIVES 

Our  objective  in  this  project  is  to  gain  increased  understanding  of  the  lateral  variations  of  velocity  structure  in  the 
crust  and  upper  mantle  of  the  region  expressing  the  continental  collision  between  the  Arabian  and  Eurasian  plates. 
Our  strategy  is  to  combine  different  types  of  geophysical  data  sets,  in  this  case,  body  wave  travel  time  data  and 
gravity  data,  in  a  formal  joint  inversion  for  velocity  structure.  Preliminary  work  has  shown  that  body  wave  data  and 
gravity  data  can  be  combined  to  yield  improved  resolution  of  the  variability  of  crustal  and  upper  mantle  velocities.  A 
key  element  to  this  success  is  careful  filtering  of  the  gravity  data  so  that  the  signal  source  region  corresponds  to  the 
depth  range  that  is  sampled  by  the  seismic  data. 

Our  choice  of  region  for  this  study  is  driven  in  part  by  the  availability  of  an  exceptional  data  set  of  body  wave  travel 
time  data  that  results  from  a  series  of  research  projects  on  calibrated  locations  in  the  region.  No  other  area  of 
monitoring  interest  is  so  well  sampled  by  ground  truth  events.  This  project  is  one  effort  to  make  maximum  use  of 
this  unique  data  set. 

RESEARCH  ACCOMPLISHED 

We  have  made  progress  in  all  major  areas  of  this  project,  acquiring  the  earthquake  and  gravity  data  sets, 
development  of  the  joint  inversion  code,  and  performing  preliminary  inversions. 

Calibrated  Earthquake  Dataset 

We  have  imported  a  large  data  set  of  calibrated  earthquake  locations  and  associated  station  information  and  arrival 
times  into  the  Ground  Truth  Database  that  is  maintained  at  Los  Alamos  National  Laboratory  (LANL).  This  data  set 
was  developed  in  several  previous  projects  in  the  Nuclear  Explosion  Monitoring  Research  and  Development 
(NEM  R&D)  research  program  (Bergman  et  al.,  2009).  These  are  all  available  calibrated  events  for  the  study  region 
at  this  time,  although  we  expect  to  be  able  to  add  some  additional  calibrated  clusters  during  the  project. 

By  “calibrated”  we  mean  that  the  hypocenters  have  been  estimated  with  a  specialized  multiple  event  relocation 
process  that  removes  systematic  location  errors  from  clusters  of  earthquakes,  by  establishing  the  location  of  the 
cluster  using  only  near-source  data  that  are  insensitive  to  unknown  velocity  structure.  Such  data  includes  seismic 
recordings  at  short  epicentral  distances,  field  mapping  of  faulting,  and  remote  sensing  analyses  of  ground 
displacements  such  as  Interferometric  Synthetic  Aperture  Radar  (InSAR).  The  accuracy  of  relative  locations  of 
clustered  events  is  also  improved  in  the  relocation  analysis.  The  majority  of  events  in  the  data  set  qualify  as  GT590 
and  many  qualify  at  lower  levels  of  location  accuracy.  In  the  relocation  process  the  arrival  time  data  sets  are 
groomed  to  remove  outlier  readings  and  to  estimate  empirical  reading  errors. 

The  import  process  required  extensive  review  of  station  codes  and  coordinates  to  maintain  the  correct  associations 
with  station  information  and  associated  arrival  time  data  already  in  the  database.  The  calibrated  earthquake  data  set 
is  now  easily  accessible  for  input  to  the  joint  inversion  code  that  is  maintained  at  LANL,  and  it  is  also  easy  to 
incorporate  other  (non-calibrated)  earthquake  data  in  the  region  in  the  inversion  as  needed.  The  calibrated 
earthquake  data  set  includes  1902  events  recorded  at  376  stations,  yielding  26,043  P  readings  and  6,390  S  readings 
that  we  have  used  in  the  joint  inversion.  The  distribution  of  calibrated  events  is  shown  in  Figure  1  and  the  ray 
coverage  for  the  Pn  phase  is  shown  in  Figure  2. 
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Figure  1.  Map  of  the  study  region,  showing  the  locations  of  earthquakes  with  calibrated  locations  that  are 
used  in  this  study.  Calibrated  relocations  have  been  done  with  a  multiple  event  analysis  applied  to 
clusters.  The  bounding  boxes  of  the  clusters  are  shown. 
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Figure  2.  Pn  raypaths  in  the  calibrated  earthquake  data  set.  The  data  set  also  contains  Sn  arrivals,  as  well  as 
crustal  (Pg,  Sg)  and  teleseismic  arrivals.  The  raypaths  shown  are  summary  rays  from  individual 
clusters  (Figure  1),  shown  as  red  circles,  to  stations  (white  triangles).  Most  raypaths  are  sampled 
many  times,  providing  useful  information  on  variability  of  the  input  data. 
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Gravity  Data 

We  used  gravity  observations  extracted  from  the  global  gravity  model  derived  from  the  Gravity  Recovery  and 
Climate  Experiment  (GRACE)  satellite  mission.  A  detailed  description  of  the  satellite  data  and  the  computation  of 
the  gravity  field  are  given  by  Tapley  et  al.  (2005).  The  new  gravity  field  is  more  accurate  than  previous  data  and 
reveals  previously  unobserved  tectonic  features.  Figure  3  shows  the  free-air  gravity  anomalies  for  the  study  region. 
Blue  colors  represent  gravity  highs  and  gravity  lows  are  represented  in  red.  Free-air  gravity  anomalies  contain 
information  not  only  of  the  subsurface  density  but  also  topography.  While  in  flat  regions  this  may  not  be  a  problem, 
in  areas  of  high  topographic  relief  this  effect  needs  to  be  removed;  thus,  we  converted  free-air  anomalies  into 
Bouguer  anomalies  assuming  a  standard  density  for  crustal  rocks  of  2670  kg/m3. 

To  increase  the  usefulness  of  the  gravity  data  for  investigating  3-D  seismic  structure,  we  are  exploring  high-pass 
gravity  filtering.  Though  commonly  applied  in  gravity  forward  modeling  (e.g.,  Simiyu  and  Keller,  1997),  such 
techniques  are  not  usually  applied  in  joint  inversion  studies  (e.g.,  Zeyen  and  Achauer  1997;  Tiberi  et  al.,  2003; 
Maceira  and  Ammon,  2009).  Filtered  gravity  anomalies  are  expected  to  possess  their  highest  resolving  power  at 
short  wavelengths  and  thus  enhance  resolution  of  lithospheric  structure.  Maceira  et  al.  (201 1)  have  shown  the 
possible  benefits  of  filtering.  After  high-pass  filtering,  the  remaining  short- wavelength  signal  can  be  assigned  with 
greater  confidence  to  structures  within  the  crust,  reducing  the  mutually  degrading  effects  of  smearing  between 
crustal  and  mantle  structures.  To  high-pass  filter  the  Bouguer  gravity  data  we  follow  the  filtering  parameter  choices 
of  Tessema  and  Antoine  (2004),  who  use  an  upward  continuation  method  and  demonstrate  correlation  with  crustal 
geology. 


-150  -100  -50  0  50  100  150 

Free-air  Gravity  Anomaly  (mgal) 


Figure  3.  Free  Air  gravity  anomaly  field  in  the  study  area. 

Joint  Inversion  Code 

During  a  previous  project  in  the  NEM  R&D  research  program  and  in  collaboration  with  MIT,  we  developed  an 
algorithm  named  JointTomoDD  (Maceira  et  al.,  2010),  which  is  a  modification  of  the  Maceira  and  Ammon  (2009) 
joint  inversion  code,  in  combination  with  the  double-difference  tomography  program  tomoFDD  (Zhang  and 
Thurber,  2003,  2006),  with  a  fast  LSQR  solver  operating  on  the  gridded  values  jointly.  The  model  representation  is  a 
combination  of  columns  of  rectangular  prisms  (needed  for  the  gravity  forward  modeling)  embedded  in  a  grid  whose 
nodes  are  located  at  the  center  of  each  prism. 
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Preliminary  Inversions 

We  first  performed  several  runs  using  only  the  seismic  body  waves  (P  and  S)  arrival  times.  The  resulting  3D  model 
depends  strongly  on  the  initial  model.  Figure  4  shows  the  P  and  S  velocity  fields  obtained  when  using  akl35  as  the 
starting  velocity  model.  Figure  5  shows  the  corresponding  Moho  isosurface  (choosing  a  velocity  of  8  km/s)  in  which 
quite  unrealistically  deep  Moho  depths  can  be  seen.  Figure  6  shows  the  velocity  fields  obtained  by  starting  with  a 
modified  akl35  model  in  which  the  Moho  is  set  to  -  50  km  depth  (as  compared  to  35  km  in  akl35),  which  is  known 
to  be  a  more  reasonable  average  Moho  depth  for  this  region.  Figure  7  shows  the  corresponding  Moho  isosurface, 
which  is  still  too  deep  in  some  places. 
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Figure  4.  Results  of  inversion  of  the  seismic  data  only  (P  velocities  on  the  left,  S  velocities  on  the  right),  shown 
at  depths  of  5, 13,  24,  52,  82  and  105  km.  The  starting  model  is  akl35. 
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Figure  5.  Moho  isosurface  (at  8  km/s)  derived  from  the  velocity  model  of  the  inversion  using  only  seismic  data 
(Figure  4). 
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Figure  6.  Results  of  inversion  of  the  seismic  data  only,  with  a  modified  version  of  akl35  (Moho  at  50  km)  as 
the  starting  model.  Otherwise,  as  Figure  4. 


Figure  7.  Isosurface  of  the  Moho  derived  from  the  inversion  of  seismic  data  alone,  using  a  modified  version  of 
akl35  (with  Moho  at  50  km)  as  the  starting  model  (Fig.  6). 

In  either  case,  the  fit  of  the  gravity  field  inferred  from  the  final  velocity  model  to  the  observed  gravity  field  is  very 
poor  (Figure  8).  In  fact  there  appears  to  be  a  significant  degree  of  anti-correlation. 

We  then  performed  some  initial  joint  inversions  combining  the  seismic  data  with  the  free-air  gravity  anomalies  for 
the  region  and  using  the  modified  version  of  akl35  as  the  initial  model.  Figure  9  shows  the  compressional  (left)  and 
shear- wave  (right)  velocity  models  from  the  joint  inversion.  The  Moho  isosurface  derived  from  this  inversion  is 
shown  in  Figure  10.  These  results  are  very  preliminary  so  we  do  not  attempt  an  interpretation.  Even  though  the  fit  to 
the  gravity  data  has  greatly  improved  (compare  Figures  8  and  1 1),  we  can  already  see  signs  of  leakage  of  the  high 
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wavenumber  gravity  signal  into  the  upper  mantle  velocity  structure.  Gravity  filtering  of  the  Bouguer  anomalies  will 
help  alleviate  this  problem  and  we  are  now  in  the  process  of  performing  those  runs. 


Figure  8.  Comparison  of  the  gravity  field  predicted  from  the  velocity  model  shown  in  Figure  6  with  the 
observed  gravity  field.  There  is  a  significant  degree  of  anti-correlation. 
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Figure  9.  Results  of  the  joint  inversion  of  seismic  and  gravity  data  (P  velocities  on  the  left,  S  velocities  on  the 
right),  shown  at  depths  of  5, 13,  24,  52,  82  and  105  km.  The  starting  model  is  akl35,  modified  to 
have  a  50  km  thick  crust. 
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Figure  10.  Isosurface  of  the  Moho  (8  km/s)  based  on  the  velocity  field  from  joint  inversion  of  seismic  and 
gravity  data  (Figure  9). 


Figure  11.  Comparison  of  the  gravity  field  predicted  from  the  velocity  model  shown  in  Figure  9 

(joint  inversion  of  seismic  and  gravity  data)  with  the  observed  gravity  field.  The  agreement  is, 
unsurprisingly,  much  better  than  the  case  where  the  velocity  model  is  determined  using  only 
seismic  data  (Figure  8). 

In  our  modeling  we  couple  seismic  and  gravity  data  through  a  density- velocity  relationship.  This  is  an  ambiguous 
relationship,  however,  because  seismic  wavespeeds  are  governed  not  only  by  rock  density  but  also  by  bulk  and  shear 
moduli,  which  vary  independently  and  are  influenced  by  many  factors,  including  lithology,  temperature  and 
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lithostatic  pressure.  Our  initial  inversions  have  used  a  combination  of  traditional  empirical  relationships:  one 
appropriate  for  sedimentary  rocks  (after  Nafe  and  Drake,  1963)  and  one  for  basement  rock  (Birch,  1961).  Both 
require  the  additional  assumption  of  a  fixed  Poisson  ratio  in  order  to  obtain  shear  wave  velocity.  In  subsequent 
modeling  we  are  employing  a  shear  velocity  vs.  density  relationship  proposed  by  Ludwig  et  al.  (1970),  as 
implemented  in  a  code  provided  by  D.G.  Harkrider.  This  method  allows  us  to  treat  Poisson’s  ratio  as  a  variable  in 
the  joint  inversion,  although  the  degree  to  which  the  data  will  permit  this  variable  to  be  resolved  needs  to  be 
explored. 

CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  set  up  the  seismic  (body  wave  travel  times)  and  gravity  data  sets,  finished  modification  and  testing  of  the 
code  for  the  joint  inversion,  and  begun  running  inversions  to  learn  how  best  to  combine  these  types  of  data  to  infer 
features  of  the  3-D  velocity  structure  of  a  continental  collision  zone.  We  have  performed  inversions  using  only  the 
seismic  data,  finding  that  they  are  inadequate  to  accurately  image  the  3-D  velocity  field  in  the  study  region,  except 
perhaps  for  a  few  areas  with  especially  good  coverage.  Adding  the  unfiltered  gravity  data  to  the  inversion  leads  to 
results  that  are  far  more  consistent  with  the  observed  gravity  field,  as  expected,  but  much  work  will  be  needed  to 
explore  issues  of  resolution  and  uniqueness  for  the  seismic  velocity  fields  determined  in  this  type  of  inversion. 

We  are  also  exploring  the  usefulness  of  filtering  the  gravity  data  to  create  data  sets  whose  signal  is  dominated  by 
certain  ranges  in  depth  that  can  be  associated  with  specific  seismic  phases  and  epicentral  distance  ranges.  Our  code 
is  also  able  to  incorporate  seismic  surface  wave  data  and  we  will  consider  this  as  a  way  to  improve  the  sampling 
coverage  by  seismic  data  in  the  study  region. 
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